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We investigate the superfluid properties of disordered double layer graphene systems using the 
non-equilibrium Green's function (NEGF) formalism. The complexity of such a structure makes it 
imperative to study the effects of lattice vacancies which will inevitably arise during fabrication. We 
present and compare room temperature performance characteristics for both ideal and disordered 
bilayer graphene systems in an effort to illustrate the behavior of a Bose-Einstein Condensate in the 
presence of lattice defects under non-equilibrium conditions. We find that lattice vacancies spread 
throughout the top layer past the coherence length have a reduced effect compared to the ideal case. 
However, vacancies concentrated near the metal contacts within the coherence length significantly 
alter the interlayer superfluid transport properties. 



I. INTRODUCTION 

After the foundations of superconductivity were laid 
more than fifty years ago, the condensation of excitons, 
bosonic molecules comprised of one electron and one 
hole, in semiconductors became a topic of intense inter- 
est to physicists 1 . In the subsequent years, the obser- 
vation of a Bose-Einstein condensate (BEC) in semicon- 
ductors has been limited by the fact that directly bound 
excitons have such a short lifetime 2 . In recent years, the 
search for BEC in semiconducting systems has found 
significant experimental 3 13 and theoretical progressSSEZl 
in coupled quantum well systems in the Quantum Hall 
regime, where the each of the layers has a filling factor 
of viayer = ^ for a total filling factor of v tota i = 1. De- 
spite significant efforts, however, there has yet to be any 
conclusive proof that this system does exhibit excitonic 
superfluidity. This is due to the fact that in semiconduc- 
tor quantum wells, a magnetic field is required to drive 
the phase transition and this magnetic field creates edge 
states which flow along the edges of the system. There- 
fore, an argument may be made which states that it is 
just as likely that the injected quasiparticles are simply 
transported in the dissipationless edge states rather than 
in the bulk of the system as one would expect for a BEC. 

Spatially separated monolayers of graphene have 
been predicted to exhibit excitonic superfluidity at tem- 
peratures approaching room temperature^ 3 -^ before 
thermal fluctuations break the individual excitons and 
destroy the condensate. The necessary conditions re- 
quired to observe excitonic superfluidity in double layer 
graphene are interlayer separations on the order of 
1 Tim and an equal number of electrons and holes in 
the respective layers. Double layers of graphene are 
suited to the observation of broken symmetry states at 
much higher temperatures as compared to that of the 
semiconductor bilayers due to the fact that it is atom- 
ically two-dimensional, thereby significantly reducing 
the screening effects, and because it exhibits a band 
structure that is both linear and gapless which implies 
particle-hole symmetry. Experimentally speaking, this 
is a much more ideal system in which to observe exci- 



tonic superfluidity as double layers of graphene may 
nest their Fermi surfaces without the need of a mag- 
netic field to quench the kinetic energy of the system. 
This eliminates the edge states which may be obstruct- 
ing the clear observation of the BEC in the semiconduc- 
tor layers. However, it is not at all clear that the con- 
clusion that broken symmetry states can exist at tem- 
peratures at and above room-temperature in double 
layer graphene is correct. Estimates for the Kosterlitz- 
Thouless temperature obtained using a mean-field the- 
ory linearized version of the critical temperature (T c ) 
combined with static Thomas-Fermi screened interlayer 
interactions have been shown to produce a much lower 
transition temperature^. This is a result of the fact that 
graphene has more fermionic degrees of freedom which 
serve to screen out the interlayer Coulomb interactions 
which will drive the phase transition. 

In situations in which such discrepancies in the loca- 
tion of the phase boundary exist between theories, ex- 
periment will be the ultimate arbiter. While the benefits 
of studying such a system experimentally are clear, the 
issues surrounding the fabrication of such an intricate 
system are plentiful. As such, it is highly unlikely that 
the graphene layers which comprise our double layer 
system will be of pristine quality. Current efforts to fab- 




FIG. 1. Schematic depiction of the system we study in this 
work. 
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ricate this system have shown tantalizing clues which 
demonstrate that the transport properties of separated 
layers are indeed modified by the interplay of inter- 
layer and intralayer Coulomb interactions 32 . However, 
the transport characteristics also do not exhibit behavior 
consistent with the formation of an exciton superfluid. 
Disorder introduced when adding the top graphene 
layer may well be the culprit which is currently occlud- 
ing the possible phase transtion as is evidenced by the 
marked reduction in mobility of the top graphene layer 
as compared to the bottom layer. While there are numer- 
ous theoretical 33 39 and experimenta l 40 ! 41 ' studies which 
outline the significant reductions in transport properties 
of single graphene layers when even small amounts of 
disorder are present in the system, there are currently 
few which consider the role of disorder in the case of 
double layer graphene. There are predictions which 
describe the role of disorder in reducing T c in double 
layer graphene 42 but none which detail how the inter- 
layer and intralayer transport in double layer graphene 
is modified in the presence of varying amounts of disor- 
der present. 

In this paper, we self-consistently solve the equa- 
tions for electrostatics in conjunction with the non- 
equilibrium Green's function formalism (NEGF) and lo- 
cal interlayer exchange interactions to understand the 
evolution of separately contacted double layer graphene 
in the regime of excitonic superfluidit y 43 " 45 ! We study 
this system under finite bias with varying amounts and 
locations of disorder. We limit the disorder in our study 
to be vacancies in the top graphene layer while the bot- 
tom layer is assumed to be ideal in order to better un- 
derstand the situation present in the experiments. We 
begin in Section |ll] by outlining the methods and ap- 
proximations we employ to illustrate interlayer quan- 
tum transport near the proposed regime of excitonic su- 
perfluidity in systems where one of the layers contains 
disorder. In Section III we present the results of our self- 
consistent quantum transport theory for an ideal system 
of nanometer sized zig-zag graphene sheets with perfect 
A-A registration between the top and bottom layers. We 
focus on the calculation and analysis of the critical tun- 
neling currents^ 11 ! 45 ! 46 !, which denote the end of coherent 
interlayer transport and are marked by a significant rise 
in the interlayer resistance for currents beyond the crit- 
ical value, in agreement with previous prediction d 44 l 45 i 
In Section IV we present the main results of our paper 
which are our numerical calculations for the disordered 
graphene system. In the calculations we present here, 
we treat the disorder as lattice vacancies present only in 
the top graphene layer. In particular, we focus on the 
behavior of superfluids in double layer graphene where 
vacancies are present in two distinct regions in our sys- 
tem relative to the coherence length, or the distance over 
which wavefunctions penetrate into the insulating su- 
perfluid state: (i) in the center of the layer past the co- 
herence length and (ii) within the coherence length of 
the contacts used to inject and extract current. When 



disorder is only included in the regions past the coher- 
ence length in the top layer, we find that the critical cur- 
rent is reduced and follows a square root dependence 
on the top layer vacancy concentrations. However, as 
the disorder only makes local perturbations to the sys- 
tem, reducing the local exciton concentrations, we find 
that the the critical current is only degraded by 30% of 
its original ideal value. In contrast, when vacancies lie 
within the coherence length, the corresponding reduc- 
tion in available area of conductive graphene limits in- 
terlayer conductance. This results in an imbalance in 
superfluid current on each side of the system, yielding a 
critical current less than 80% of its ideal value and which 
depends linearly on the concentration of vacancies in 
the top layer. Furthermore, for both disorder cases we 
consider here, we find that no steady-state superfluid 
density can be found when more than 4% vacancies are 
included in one monolayer. 



II. SIMULATION METHODOLOGY 

In Fig. [TJ we plot a schematic of the system of interest. 
In this work, we consider two zigzag graphene mono- 
layers assumed to be perfectly aligned with one another 
and separated by a thin dielectric. Contacts along the 
edges of each of the layers inject and extract current, 
and top and bottom gates manipulate the quasiparti- 
cle concentrations in each of the layers. In this way, we 
may tune them to contain the proper quasiparticle con- 
centrations predicted by many-body theory 2 — 30 to in- 
duce a superfluid phase transition. The top and bottom 
gates are separated from the graphene layers by 20 nm 
of SiGv The two graphene layers are separated from 
one another by 1 nm SiC>2 spacer dielectric to be in the 
regime of superfluidity predicted by many-body calcu- 
lations. We consider the oxide regions to be perfect in 
the sense that they do not contain any stray charges and 
have perfect interfaces with the graphene layers. We 
choose the x direction to lie along the length of the sys- 
tem, the y direction to lie along the width and the z along 
the depth. We choose each monolayer to be 30 nm long 
by 10 nm wide. We choose these dimensions rather than 
larger system sizes so that we can calculate the atomistic 
transport properties of an increased number of disorder 
distributions within a reasonable amount of time. The 
top and bottom gates (Vtg = —Vbg) are are gated to 
effect individual carrier concentrations of 10 13 cm~ 2 in 
each layer, which corresponds to a Fermi energy in the 
top (bottom) layer of 0.4eV (-0.4eV). The gate bias con- 
ditions and the interlayer separation has been chosen so 
as to satisfy the conditions for room temperature pseu- 
dospin ferromagnetism 28 . Here, we focus on the regime 
where the layer electron and hole populations place the 
system firmly in the dense electron-hole regime where 
we expect the electron-hole pairs to form a BCS-type 
state. This is as opposed to the dilute limit of electron- 
hole densities which can be described as a weakly inter- 
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acting Bose system of excitons®. As we wish to collect 
a sufficient statistical distribution on the effects of va- 
cancy distribution on the interlayer transport properties 
of double layer graphene, we utilize several randomized 
vacancy configurations for each concentration we exam- 
ine in this work. 

We begin the NEGF simulatio n 1 48 ! 49 ! with the atom- 
istic tight-binding description of an individual graphene 
monolayer, 



Htl = Y, r\t)(j\+V l \ l )(t\, 



(1) 



<hJ> 



where lattice points i and j are first nearest neighbors. 
t = — 3.03eV is the nearest neighbor hopping energy for 
the p z orbital of graphene, which allows for the unique 
low-energy linear dispersion at the K and K' points in 
the Brillouin zone. We neglect hopping among further 
nearest neighbors and other orbitals, as nearest neigh- 
bor p z orbital hopping is the predominant interaction for 
graphene in the probed energy range. The on-site poten- 
tial energy Vi = (f>{vi) is calculated via a 3-dimensional 
Poisson solver. We use a phenomenological model to 
simulate a generic metal contact with a constant density 
of state d 50 ! 51 ! . This model captures the basic self-energy 
needed to appropriately simulate a metal contact with- 
out taking into account multiple orbitals or complex in- 
terface problems such as lattice mismatch or Schottky 
barrier height. 

We may now generalize our layer Hamiltonian to 
the double layer Hamiltonian by coupling the top and 
bottom monolayers with the following Bogoliubov-de 
Gennes (BdG) Hamiltonian, 



T~iBdG 



H TL 
H BL 



(2) 



with the interlayer interactions including both single 
particle tunneling and the mean-field many-body con- 
tribution, A, coupling the two layers using a local den- 
sity approximation. In Eq. j2|, \i represents a vector that 
isolates each of the Cartesian components of the pairing 
vector, a M represents the Pauli spin matrices in each of 
the three spatial directions, and ® represents the Kro- 
necker product. 

In order to correctly account for the dynamics of 
the double layer graphene system, we first include the 
many-body interlayer interactions in the Hamiltonian. 
Within the Hartree-Fock mean-field approximation, we 
may define interlayer interactions through the expecta- 
tion value of the full Hamiltonian, 



(til H BdG 



(3) 



We assume that the graphene monolayers are perfectly 
registered, so that electrons at site i in the top layer ((f» | ) 
only bind with holes at site j on the bottom layer 
when i = j. U is the strength of the interlayer on-site 
Coulomb interaction, whose selected value we address 
later. 



m exc is the magnitude of the order parameter result- 
ing from our analysis of BECs of indirectly bound ex- 
citons. It is proportional to the off -diagonal terms in 
t he sin gle particle density matrix, which is represented 

a j48l52l 
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The on-diagonal density matrix (pff, pn) corresponds to 
the associated electron and hole densities of the top and 
bottom monolayers, respectively. The order parameter 
m exc can now be defined as a function of the interlayer 
component of the density matrix, 



< c = Pn + Pit = 2Re (Pn)> 

m exc = -iPn + iPll = 2Im{Pn)- 



(5) 



The density matrix is directly calculated within the 
NEGF formalism. It is consequently both an input and 
an output of our simulation. We iterate over the above 
mean-field equations, in conjunction with the Poisson 
equation for electrostatics, to obtain a self -consistent so- 
lution with compatible particle densities and potential 
profile utilizing the Broyden method 53 to accelerate con- 
vergence. 

We can now expand Eq. Q to show a simplified 
form for our BdG Hamiltonian in which interlayer in- 
teractions are expressed in terms of their directional 
components®, 



T-Lbclg = 



H 



TL 



i-A z A*- 
iA y H B l 



iA, 



(6) 



The directional components of the interlayer interac- 
tions A are expressed as 



A, = (A sas + Um x exc ) 
A y = Um v exc 



A* = jj(Vf | | -Vi | |). 



(7) 



The on-diagonal term in the interlayer interactions, A z , 
is due to screening caused by the unbound carriers in 
each monolayer, and acts to separate the two Fermi sur- 
faces. The value of the single particle tunneling energy, 
A sas , is proportional to the probability of a single elec- 
tron tunneling events through the thin dielectric and re- 
combining with a hole. Single-particle tunneling is an 
adverse event, thus it is desirable to have a very thick 
barrier in between the two graphene layers to maximize 
the lifetime of the indirectly bound excitons. However, 
we also need strong interlayer Coulomb interactions in 
order to drive the superfluid phase transition which ne- 
cessitates a thin barrier. In this paper, we set A sas = 
1/ieV, sufficiently small so that the lifetime of the indi- 
rectly bound exciton is long enough to observe conden- 
sation but not so small as to require an intractably large 
number of iterations before self-consistency is reached. 
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FIG. 2. Interlayer and intralayer current plots for the ideal sys- 
tem as a function of counterflow bias (solid black line and red 
dashed line, respectively). 



Results of previous many-body calculations show the 
value of the order parameter is approximately one tenth 
the Fermi energy 2 ^. Our simulations show that this 
value of m exc corresponds to an interlayer coupling 
strength of 2.0eV. Although this is less than the un- 
screened mean-field interaction strength 52 (U = % f=s 
4.6eV), it is more plausible as it factors in damping ef- 
fects due to screening^. We set U = 2.0eV in all simu- 
lations for the duration of this paper. We calculate the 
magnitude of the order parameter using the expectation 
values of the density matrix in Eq. |5|, 



(8) 



Similarly, we may identify the phase of the order pa- 
rameter dependent on the same expectation values of 
the density matrix 4 ^, 



tan 



(9) 



As we expect that the interlayer transport properties 
are similar to the case of a Josephson junction, we ex- 
pect 4> e.xr = 7r/2 at zero bias to maximize interlayer 



curren t The quasiparticle and condensate current 
densities, the artifacts of the Andreev reflection action, 
are proportional to the spatial phase gradient and mag- 
nitude of the order parameter^. 



III. INTERLAYER TRANSPORT IN IDEAL DOUBLE 
LAYER GRAPHENE 

It is vital to first understand the interlayer transport 
properties in ideal double layer graphene systems com- 
prised of zigzag graphene nanoribbons. After estab- 
lishing ideal system properties, we can compare trans- 
port in the disordered scenarios to see how disorder im- 
pedes the inter and intralayer current flows. To obtain a 




FIG. 3. The dispersion relation for an ideal bilayer graphene 
system in the transport direction (a) at a counterflow bias of 
Vtl — —Vtr — 0.05 V which is below the critical current and 
(b) at a counterflow bias of Vtl = — Vtr = 0.23 V which is 
above the critical current. 



more in-depth understanding of the interlayer transport 
in double layer graphene systems, we draw a compari- 
son to the Andreev reflections 54 that occurs at a metal- 
superconductor interface to explain the non-equilibrium 
physics of the condensate. In a superconducting sys- 
tem, when an electron with energy less than the super- 
conducting band gap is injected into a superconductor, 
the injected electron penetrates a certain distance before 
producing a Cooper pair that moves across the super- 
conductor and a retro-reflected hole of opposite spin in 
the normal metal. For the exciton superfluid, an elec- 
tron with energy less than the superfluid band gap in- 
jected into the system will penetrate a short distance 
in the top layer before driving an exciton that moves 
across the channel. Put in the language of our dou- 
ble layer graphene system, this process corresponds to 
an electron being injected into the top left contact that 
causes an exciton to move across the system while an 
electron is reflected into the bottom left contact to con- 
serve current 55 . 

The drag-counterflow geometry (Vtl = —Vtr;Vbl = 
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Vbr — V) will cause an electron to be reflected into 
the bottom left contact and a hole into the bottom right 
contact 52 , inducing an effective current through the bot- 
tom layer. All non-equilibrium configurations described 
below are in this geometry, with the left and right con- 
tacts on the top monolayer set to magnitude of the coun- 
terflow bias parameter. The entire process results in a 
condensate current, due to the propagating exciton, and 
a quasiparticle current, caused by the injected and retro- 
reflected individual carriers. The quasiparticle current 
is only nonzero within the coherence length (L c ), the 
maximum length an injected particle penetrates into the 
superfluid gap before triggering the exciton scattering 
event. The microscopic attributes of the superfluid and 
the dynamics of the Andreev reflection may not be read- 
ily apparent in experiments, but will have a macroscopic 
effect on the observable interlayer current. 

In Fig. |2j we present these observable interlayer and 
intralayer transport properties for an ideal system. The 
currents are odd functions of the bias, reflecting the sys- 
tem's ambipolar nature. At low bias, interlayer current 
is linear and interlayer conductance is constant, as in- 
jected carriers see the superfluid gap and trigger the 
Andreev event discussed above. Low-energy injected 
carriers are unable to pass through the superfluid gap, 
and intralayer current is negligible within this range. 
States above the superfluid band gap do begin to form 
at > 0.14 eV, as seen in the superfluid band struc- 
ture in Fig. [3^. Interlayer transmission quickly van- 
ishes when higher-energy states become accessible in 
the spectrum; transport along the monolayer dominates 
in these ranges. This is evidenced by the emergence of 
intralayer current and a drop in interlayer conductance 
beginning at a counterflow bias of 0.14 V. 

Beyond the critical counterflow bias, 0.22±0.01 V in 
the ideal system, the superfluid can no longer adjust 
its phase to accommodate the current flow. When self- 
consistency is lost, only single particle tunneling con- 
tributes to interlayer tunneling. The small magnitude 
of the resultant interlayer current, less than 1 pA for a 
bias of Vtl = —Vtr = 0.25 V, is negligible compared 
to the many-body contribution to the interlayer current 
observed when the system obtains self -consistency. Ad- 
ditionally, past the critical current, the transport charac- 
teristics of the bilayer becomes time-dependenPS and 
the mean field equations no longer represent a valid de- 
scription of the physics. Although a transient remnant 
of the superfluid may yet persist beyond the critical bias, 
the many-body portion of the Hamiltonian is approxi- 
mately zero and the time dependent oscillation in the 
interlayer transport relationship can not be captured in 
our steady-state simulations. 

The phase transition is observed in Fig. |3b, which 
shows the dispersion for the ideal system after critical 
current has been passed. There is no appreciable change 
to the band structure within each phase. When no su- 
perfluid exists, each monolayer has a vanishing band 
gap and linear dispersion. As expected, the Dirac points 



occur in the normal phase at the K' and K points in the 
Brillouin zone with a Fermi energy of ±0.4eV. The band 
gap vanishes at biases larger than this critical transition 
value. The closing of the band gap after the phase tran- 
sition allows for more low energy states, causing a spike 
in intralayer current well beyond the range of interlayer 
transport in the superfluid phase. 

The band structure in the transport direction of the 
system is derived 56 from a small portion of the con- 
verged Hamiltonian by assuming m exc is periodic in the 
transport direction. This holds true only for the ideal 
system, where m exc is smooth and consistent through- 
out the channel with a magnitude near 10% of the Fermi 
energy and a phase of tt/2. Randomly-placed impuri- 
ties and vacancies, however, break translational symme- 
try so that the Hamiltonian is no longer periodic in the 
transport direction. We are thus unable to calculate the 
dispersion relations of disordered systems. 

We find that the ideal system's critical current is I c = 
±18.4 uA, at a counterflow bias of ± 0.22 V. We may 
compare this with the analytic approximation for the 
critical current in a system where coherence length is 
smaller than system length 45 . In this case the conden- 
sate must satisfy an elliptic sine-Gordon equation, 

A 2 V 2 - sin{(f>) = 0. (10) 

When this equation is solved in the static case, we obtain 
a relatively simple expression for the critical current in 
our system, 



Therefore, for a system of width W = 10 nm with a co- 
herence length of L c w 5nm and order parameter mag- 
nitude of p s m 0.04eV, the analytic critical current is 
roughly if w 19.5/zA This is in good agreement with 
the critical current calculated in our simulations of an 
ideal bilayer. The reduction in the critical current is ex- 
pected in our system as reflections off the superfluid 
gap lead to enhanced differences in local interlayer in- 
teraction terms which enter into the Hamiltonian self- 
consistently through the A z term and serve to energeti- 
cally separate the layers! 4 ^ 

Fig. fi] plots the interlayer quasiparticle current den- 
sities along the transport direction averaged along the 
width of the system in an ideal system with a coun- 
terflow bias of 0.05V and 0.20V, just before the phase 
transition to a normal state. Quasiparticle current tun- 
neling, which is part of the same process that launches 
the exciton through the channel, is highest near the con- 
tacts where carriers are injected, and evanescently de- 
cays into the channel, so that carriers which provoke in- 
terlayer transmission are only present a small distance 
into the channel. We find this distance, called the coher- 
ence length (L c ), to be roughly 5 nm by inspection of the 
plot. 

Note in Fig. B] that interlayer quasiparticle current 
magnitude is nearly equivalent at low bias but becomes 
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FIG. 4. The interlayer quasiparticle current density along the 
transport direction in the ideal bilayer system averaged along 
the width of the system for counterflow biases of 0.05V (red 
dashed line) and 0.20V (blue solid line). The current density 
is only nonzero within the coherence length, which by inspec- 
tion is approximately 5 nm. Quasiparticle current density is 
approximately odd about the center of the system at low bias, 
but symmetry is lost when a high contact bias limits carrier 
concentration on one side of the system. 

very asymmetric near the critical transition point. The 
Poisson equation in the counterflow bias configuration 
restricts the magnitude of the condensate order param- 
eter because the negatively biased top right contact lo- 
cally decreases the density of electrons able to pair into 
the condensate. Conversely applying a bias across the 
hole-doped bottom layer will limit superfluid density 
and quasiparticle current near the left (positively biased) 
contact. The decreased density results in a smaller in- 
terlayer current density on one side, as seen in the case 
with higher counterflow bias. All interlayer currents 
are the currents generated by the contacts on the right 
side of the system, as this is the side with limited exci- 
ton density is the transport bottleneck. When the dis- 
parity in condensate currents on each side of the bilayer 
reaches the critical limit, steady-state supercurrents are 
no longer possible^. 

IV. INTERLAYER TRANSPORT IN DISORDERED 
DOUBLE LAYER GRAPHENE 

We now seek to understand how the behavior of the 
interlayer transport properties change as we introduce 
disorder into the system. We show the main result of 
our analysis in Fig. |5j which plots the statistical varia- 
tion in critical current as a function of vacancy concen- 
tration. In this figure, we show two distinct types of be- 
havior characteristic of the two categories of locations of 
top layer vacancies we consider. In the middle of the 
figure, at a vacancy concentration of 0%, we plot the 
critical current we obtain for the ideal case. To the left 
of the ideal case, we plot the statistical variation in the 



FIG. 5. Statistical variation of critical interlayer current for sev- 
eral runs at the bias just before the condensate is broken. Va- 
cancies which are more than a distance of one coherence length 
away from the top layer contacts are to the left of the ideal 
(0%) case, depicted by a black circle, while vacancies which 
are within one coherence length from the top layer contacts 
are to the right. The excitonic condensate is lost in both cases 
for top layer vacancy concentrations larger than 4%. The solid 
lines represent analytic calculations using Eq. (TT) with modi- 
fied values obtained from calculations as described in Section 

El 

critical current for random top layer vacancy concentra- 
tions which are beyond the coherence length, L c , from 
either contact on the top layer. In this case, we see that 
as we increase the random top layer vacancy concentra- 
tions, there is a decrease in the critical current but that 
the decrease is proportional to the square root of the top 
layer vacancy concentration. Furthermore, we also no- 
tice that the mean of the statistical distribution of critical 
currents associated with various top layer vacancy dis- 
tributions varies only by approximately 5% even as the 
concentration of vacancies is increased. 

As we move away from the ideal case towards the 
right hand side of Fig. |5j we plot the statistical change 
of the critical current associated with within a distance 
of L c from either the left or right contact of the top layer. 
Here we see quite different behavior as compared to the 
situation where the vacancies are more than L c from the 
top layer contacts. Now the critical current has a very 
distinct linear decrease as we increase the vacancy con- 
centrations. We also see that there are very large vari- 
ations in the location of the critical current as we shift 
the locations of the top layer vacancies. In the subse- 
quent sections, we will explore the physics of these two 
distinct regions and explain the observed dependence of 
the critical current seen in each situation. 



A. Layer Disorder beyond L c 

We model disorder as randomly placed vacancies 
within the top layer of our double layer system. To 
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Counterflow Bias (V) 

FIG. 6. Set of I-V curves for systems with varying concentra- 
tions of randomly placed vacancies at least 5 nm away from 
the top layer contacts. 



model a vacancy, we modify the ideal Hamiltonian by 
setting all hopping to a missing atom to zero. This ef- 
fectively blocks any interaction with the vacancy by set- 
ting the tight-binding overlap integral of the spatial p z 
orbital states to zero^A We randomly remove a fixed 
percentage of carbon atoms from the specified region of 
the top monolayer, leaving the bottom monolayer un- 
perturbed, and perform our numerical calculation. 

In Fig. [6] we plot the interlayer current as a function 
of voltage oias for various concentrations of vacancies in 
the top layer of our double layer graphene system along 
with the ideal double layer case for comparison. We find 
that there is little reduction in the interlayer conductiv- 
ity regardless of vacancy concentration or vacancy loca- 
tion past the coherence length at low counterflow bias. 
The phase transition to two normal, incoherent Fermi 
liquids, however, occurs at a lower bias than in the ideal 
case. Channel disorder decreases critical current by 20% 
for 1% vacancy concentration to 30% for vacancy con- 
centrations of 4%. Beyond top layer vacancy concentra- 
tions of 4%, we find no self -consistent solution and the 
interlayer transport is dominated by single particle tun- 
neling events. Resultant interlayer currents are orders of 
magnitude smaller as the majority of the current injected 
into the system now flows across the graphene layers. 

In order to better understand these reductions in in- 
terlayer current as the top layer vacancy concentration 
is increased, we will examine the magnitude of the or- 
der parameter, \m exc \. In Fig. [7j we plot the magnitude 
of the order parameter for a random vacancy concentra- 
tion of 1% . In this situation, the order parameter magni- 
tude does not remain constant over the entire system at 
10% Ef, as was the case in the ideal system. The vacan- 
cies locally destroy the condensate 31 and reduces \m exc \ 
at surrounding points up to a distance of 0.5 nm, or three 
to four nearest neighbors from the vacancies. However, 
we find no long range effect is seen when vacancies are 
isolated from one another by more than approximately 
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FIG. 7. Magnitude of the order parameter for a system with 
1% vacancies in the channel. The condensate is very close to 
that of the ideal case, with a drop in magnitude only occuring 
locally near disorder. As expected, |m e3; c| ~ O.lEp away from 
disorder. 




Length (nm) 

FIG. 8. Magnitude of the order parameter for a system with 4% 
vacancies in the channel. Magnitude clearly drops nonlinearly 
with the increase in disorder, but quasiparticle current is still 
able to form near the contacts. 



2 nm. At this particular vacancy concentration, we find 
that m exc is reduced by 40% over an appreciable area of 
our system as a result of these vacancies. 

This situation is to be contrasted with Fig. [8] where 
we plot m exc for a random top layer vacancy concentra- 
tion of 4%. In this case, significant areas clearly emerge 
where concentrated disorder has long range effects on 
superfluid density. These areas of high vacancy concen- 
trations in Fig. 8^ give rise to values for |m e2;c | that are 
less than 20% the ideal value over significant areas of 
the system. Superfluid magnitude remains at its ideal 
level near the contacts, sufficiently quarantined from the 
vacancies. 

Nevertheless, the root cause of the sublinear behav- 
ior that we see in the critical current as we increase the 
top layer vacancy concentration is not yet resolved, from 
Eq. dill. To explain the sublinear behavior, we exam- 
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ine the localized density of states (LDOS^Pl Vacancies 
induce a LDOS similar to the case of an impurity^ in 
graphene layers. In Figs. |9| and [10| we plot the LDOS 
for the ideal case and for the case of 4% top layer va- 
cancy concentration. Clearly we can see a stark contrast 
in the low-energy LDOS in the superfluid phase when 
vacancies are introduced. The LDOS closely resembles 
those of monolayer graphene at higher energies^, with 
peaks at E = Ep ± r. The ideal top and bottom mono- 
layers exhibit perfectly antisymmetric densities of states 
so that equivalent carrier concentrations arise in the op- 
positely gated top and bottom monolayers. We find that 
no states exist in the superfluid gap in the ideal bilayer, 
as one would expect for a condensate in which all of the 
quasiparticles participate. 

Disorder changes the transport properties of the sys- 
tem by introducing mid-gap states in the top layer seen 
in Fig. 10 Localized states also arise, to a lesser degree, 



in the bottom layer due to charge pileup induced by 
Coulomb attraction through the thin spacer dielectric. 
The biased top and bottom gates, necessary to generate 
the sufficient carrier concentrations, cause the undesir- 
able occupation of the mid-gap states up to the Fermi en- 
ergy. As the electrons and holes make their way across 
the layers, they scatter off of these localized states which 
changes the interlayer phase relationship between the 
top and bottom layers, which manifests itself as changes 
in the A y term in Eq. (|5j. As a result, when the decreased 
interlayer phase relationship is input into the calcula- 
tion of m exc in Eq. |8|, m exc drops from its ideal value 
of 0.041 eV to 0.029 ± 0.002 eV at 4% top layer vacancy 
concentration, a decrease of roughly 30%. The result of 
this scattering is then that the reduced interlayer phase 
component then enters into Eq. (IT) in m exc it does so 
through the square root thereby giving rise the square 
root dependence that we see in on the left hand side of 
Fig. [5] When the average value of the change in the in- 
terlayer phase is added into the calculation of the critical 
current in Eq. JTlJ and plotted along with our numeri- 
cal results on theleft hand side of Fig. |5j we find good 
agreement between the two values. 



B. Layer Disorder within L c 




-0.5 0.5 

Energy (eV) 



FIG. 9. Local density of states at low energy for an ideal sys- 
tem. The solid (dashed) line represents the top (bottom) layer. 
No localized states exist below a magnitude of 0.14eV for the 
superfluid. 




Energy (eV) 

FIG. 10. Local density of states at low energy for a system with 
4% vacancies in the top layer. The solid (dashed) line repre- 
sents the top (bottom) layer. Low energy states due to the 
vacancies introduces a weak but apparent LDOS in the ideal 
bottom layer as well, due to the thinness of the spacer dielec- 
tric. 



While the double layer graphene system is moder- 
ately robust to vacancies deep in the channel, this is not 
the case when the top layer vacancies occur within one 
L c of the contacts on the top layer. In Fig. 11 we plot 



the interlayer current as a function of the random top 
layer vacancy concentration where the vacancies occur 
within within one L c of the contacts on the top layer. 
We see that vacancies near the contacts can cause a sig- 
nificant drop in both conductivity and the voltage bias 
at which self-consistency is lost. Quasiparticle current 
density is now dependent on the location of disorder 
within the coherence length, as magnitude evanescently 
decays from the contact interface. Interlayer current is 



thus sensitive to the specific location and amount of dis- 
order within the coherence length. 

Fig. [12] shows the magnitude of the order parameter, 
m exc , for a case of 1% vacancies randomly distributed 
near the contacts. Quasiparticle current density drops 
with the presence of contact disorder, as it is propor- 
tional to the magnitude of the order parameter. Clearly, 
groups of vacancies appear on the left and right side of 
the contact that begin to show a nonlinear, long-range 
effect on superfluid density. Top layer vacancies clos- 
est to the contact, where quasiparticle tunneling magni- 
tude is largest, cause the biggest detriment to the magni- 
tude of interlayer current. Despite the vacancies reduc- 
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Counterflow Bias (V) Length (nm) 



FIG. 11. I-V curve for systems with varying concentrations of 
randomly placed vacancies within 5nm of the top layer con- 
tacts. 
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FIG. 12. Magnitude of the order parameter, in eV, for a system 
subjected to 1% vacancies near the contacts. A single defect 
does not affect the condensate appreciably, but several nearby 
can drop magnitude significantly. 
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FIG. 13. Magnitude of the order parameter, in eV, for a system 
subjected to 4% vacancies near the contacts. The contacts are 
too saturated with a large amount of tunneling to occur, and 
interlay er conductivity suffers. 



FIG. 14. Quasiparticle current density for a system with 4% 
contact disorder with a counterflow bias of 0.04V, just before 
the phase transition. No condensate exists near the contacts, so 
the largely degraded quasiparticle tunneling that does occur 
only happens beyond the disorder. 



ing the space in which quasiparticles may be injected 
without scattering, we still find that the condensate is 
able to form, and that sites within three nearest neigh- 
bors of a single defect are not appreciably affected. We 
find that little randomized bunching occurs in the 1% 
case, as the disorder is too sparse to create significantly 
different scenarios. Interlayer current remains rather ro- 
bust, roughly 30% smaller than ideal. The random na- 
ture of the placement, however, causes a high variance 
in interlayer current in this scenario which is discussed 
below. 

Larger deviations from the ideal interlayer critical cur- 
rent are found as we move to higher top layer vacancy 
concentrations. We find significant amounts of cluster- 
ing occur which leads to a significant reduction in the 
available space for quasiparticles may be injected, as 
seen in Fig. [l3j and generates a very different transport 
relationship compared to the case of channel disorder. 
Carriers injected into the system see a significant reduc- 
tion in the area in which quasiparticle tunneling can oc- 
cur near the contact, significantly reducing conductiv- 
ity. As we know that there is a linear dependence on the 
width of the system in Eq. (IT) , this gives a very simple 
explanation for the physics or the linear decrease in the 
critical current observed in Fig. [5] This information al- 
lows us to conclude that graphene with vacancies within 
the coherence length effectively reduces the width over 
which quasiparticle tunneling occurs, and causes a lin- 
ear decay in I c with respect to contact disorder strength 
until the condensate vanishes. 



In Fig. 14 we see how the average interlayer quasi- 
particle current density qualitatively shifts when top 
layer vacancies are included within the coherence length 
of the top layer contacts. Equivalent disorder concen- 
trations generate disparate condensate currents on each 
side because vacancies are randomly configured to be 
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more closely lumped near one contact than the other. 
Whereas superfluid excitons are able to permeate disor- 
der as long as the condensate exists, an increasing pro- 
portion of bare electrons and holes cannot penetrate the 
contact disorder to reach the condensate. The critical 
bias at which the phase transition occurs decreases as 
disorder increases because of the discrepancy in conden- 
sate current from each contact. The on-site potential is 
greater at vacancies near the contacts due to Vtr and 
Vt l ■ The LDOS is thus more quickly occupied, and fur- 
ther accelerates the phase transition. As a result, critical 
current roughly drops by an average of 30%, 45%, 60%, 
and 80% for vacancy concentrations of 1%, 2%, 3%, and 
4%, respectively. We can see the accuracy of our conclu- 
sion on the right hand side of Fig. [5] where we plot our 
data superimposed with a model ofthe reduction of the 
critical current expected when each vacancy is assumed 
to remove 0.5 nm from the effective width of the system 
in the effected region for a given vacancy concentration. 
Here we see good agreement again between our theory 
and numerical calculations; the model falls within the 
error bars of the simulations, where deviation is the re- 
sult of non-linear randomized bunching of disorder. 

It is clear in Fig. [5] that a large variance in critical 
current exists when vacancies lie within the coherence 
length of the contacts in the top layer. This is due to 
the fact that the quasiparticle current flow, in this sys- 
tem, is concentrated close to the contacts in the middle 
of the channel, and evanescently decays into the system. 
This is best seen in Fig. 12 where we see that the mag- 



nitude of the order parameter is significantly reduced 
from the bulk value near the edge of the system both 
near and away from the contacts. Moreover, vacancies 
that are closely lumped together exacerbate the scatter- 
ing to an even greater degree because vacancies in close 
proximity to one another have a significant non-local ef- 
fect on the perturbed m exc , as seen in Fig. 13 Also 
in Fig. 13 we see the variance in critical current de- 
cays at 4% vacancies because quasiparticle tunneling is 
almost completely prevented at such a high concentra- 



tion of vacancies. As in the channel-disordered case, no 
self-consistent solution was found for vacancy concen- 
trations greater than 4%. 



V. SUMMARY AND CONCLUSIONS 

We perform self-consistent calculations to understand 
how vacancies in individual graphene monolayers can 
hinder the performance of a double layer graphene sys- 
tem in which a room temperature exciton condensate is 
predicted to form. We find that vacancies within a co- 
herence length of the contacts significantly obstruct per- 
formance by effectively reducing the width over which 
interlayer transport occurs. For a selected value of 4% 
vacancies in this region, tunneling current at a selected 
bias can drop by more than 80% compared to the ideal 
scenario. We also find that the reduced width of the sys- 
tem caused by the presence of the top layer vacancies 
produces a linear dependence on the critical current as 
the vacancy concentration is increased. Vacancies out- 
side of the coherence length have little effect on the in- 
terlayer conductivity showing a square root dependence 
of the cricital current as the vacancy concentration in 
the top layer is increased. Critical current degrades up 
to 30% due to a phase transition at smaller bias, as I c 
and we find that the reduction is due to scattering of 
layer quasiparticles from localized mid-gap states which 
modifies the average interlayer phase relationship be- 
tween the two layers. Concentrations of vacancies larger 
than 4% in one of the layers prevents the condensate 
from forming in a steady state. 
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